Velocity Structure of Self-Similar Spherically Collapsed Halos 
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Using a generalized self-similar secondary infall model, which accounts for tidal torques acting on 
the halo, we analyze the velocity profiles of halos in order to gain intuition for N-body simulation re- 
sults. We analytically calculate the asymptotic behavior of the internal radial and tangential kinetic 
energy profiles in different radial regimes. We then numerically compute the velocity anisotropy 
and pseudo-phase-space density profiles and compare them to recent N-body simulations. For cos- 
mological initial conditions, we find both numerically and analytically that the anisotropy profile 
asymptotes at small radii to a constant set by model parameters, ft rises on intermediate scales as 
the velocity dispersion becomes more radially dominated and then drops off at radii larger than the 
virial radius where the radial velocity dispersion vanishes in our model. The pseudo-phase-space 
density is universal on intermediate and large scales. However, its asymptotic slope on small scales 
depends on the halo mass and on how mass shells are torqued after turnaround. The results largely 
confirm N-body simulations but show some differences that are likely due to our assumption of a 
one-dimensional phase space manifold. 



I. INTRODUCTION 

Recent N-body simulations have revealed a wealth of 
information about the velocity structure of halos [U-Ojj]- 
However, simulations have finite dynamic range. More- 
over, it is difficult to draw understanding from their anal- 
ysis, and computational resources limit the smallest re- 
solvable radius, since probing smaller scales require using 
more particles and smaller time steps. Hence, it seems 
natural to use analytic techniques, which do not suffer 
from resolution limits, to analyze the velocity distribu- 
tions of halos. 

Numerous authors have analytically investigated the 
density profiles of halos. Work began with Gunn and 
Gott where they analyzed the continuous accretion of 
mass shells onto an initial overdensity 041] ■ This process 
is known as secondary infall. By imposing that the mass 
accretion is self-similar, Fillmore and Goldreich Q and 
Bertschinger [8(, assuming purely radial orbits, were able 
to analytically calculate the asymptotic slope of the den- 
sity profile. Since then, there have been numerous exten- 
sions, some which do not assume self-similarity, that take 
into account non-radial motions Those works that 

do not impose self-similarity can only infer information 
about the velocity dispersion using the virial theorem. 
Hence they cannot predict a halo's velocity anisotropy. 
Those works that do impose self-similarity focus only on 
the asymptotic slopes of density profiles. 

In this paper, we analytically and numerically analyze 
the velocity structure of halos using an extended self- 
similar secondary infall model [l8j . The work which in- 
troduced this extended infall model will hereafter be re- 
ferred to as Paper I. We then compare the predictions 
of our halo model to simulation results, focusing on the 
velocity anisotropy [19j and pseudo-phase-space density 

profiles ele in. 
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Density profiles do not uniquely determine a self grav- 
itating system. In order to more fully characterize dark 
matter halos, one needs to probe their phase space distri- 
butions. The velocity anisotropy and pseudo-phase-space 
density profiles are thereby useful since they complement 
density profiles by revealing additional information about 
the phase space structure of the halo. 

In section |TT] we summarize our generalized secondary 
infall model and discuss how to numerically calculate the 
radial and tangential velocity dispersions in the halo. In 
section Hill we analytically calculate the asymptotic be- 
havior of the radial and tangential kinetic energy profiles 
on small and intermediate scales. In section HVl we com- 
pare our numerically calculated anisotropy and pseudo- 
phase-space density profiles to recent N-body results and 
conclude in section IVl 



II. SELF-SIMILAR MODEL 

Here we first summarize the self-similar halo formation 
model developed in Paper I. The model is characterized 
by four parameters {n,p, B,-cu} which are described be- 
low. 

In this model, the universe is initially composed of a 
linear spherically symmetric density perturbation with 
mass shells that move approximately with the hubble 
flow. Because of the central overdensity, mass shells 
eventually stop their radially outward motion and turn 
around. The radius at which a mass shell first turns 
around, or its first apocenter, is known as the turnaround 
radius. Since the average density is a decreasing function 
from the central overdensity, mass shells initially farther 
away will turnaround later. The halo grows by contin- 
uously accreting mass shells. Mass shells are labeled by 
their turnaround time or their turnaround radius . 

Model parameter n characterizes how quickly the ini- 
tial linear density field falls off with radius (S oc r~ n ). It 
is related to the effective primordial power spectral index 
n e g (d\nP/d\nk) through n = n e g + 3 [21]. Since n c s 
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depends on scale, n is set by the halo mass. As in Paper 
I, we restrict our attention to < n < 3 so that the ini- 
tial density field decreases with radius while the excess 
mass increases with radius. Note however that n > 1.4 
corresponds to objects larger than galaxy clusters today. 
Model parameter n also sets the growth of the current 
turnaround radius: r ta oc where /3 = 2(1 + n)/3n fl8| . 

Self-similarity imposes that at time i, the angular mo- 
mentum per unit mass L of a particle in a shell at r, and 
the density p and mass M of the halo have the following 
functional forms flal. 



L(r,t) 
p{r,t) 



B- 
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M(r,t) = f PB (t)rl(t)M(X) 
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(3) 



where A = r/i\ a (t) is the radius scaled to the current 
turnaround radius, ps — l/6irGt 2 is the background den- 
sity for an Einstein de-Sitter (fiat Sl m = 1) universe, and 
B is a constant. Inspired by tidal-torque theory and nu- 
merical simulations, we take / to be: 




if t < U 
if t > t t 



(4) 



Model parameter p, defined above, sets how quickly 
angular momentum builds up before turnaround while B 
sets the amplitude of angular momentum at turnaround. 
In Paper I, using cosmological linear perturbation theory, 
we constrained p and B so that the angular momentum 
of particles before turnaround evolves as tidal torque the- 
ory predicts [22-25]. Conveniently both p and B are set 
by the halo mass. However, after comparing to density 
profiles from N-body simulations, we found that our ex- 
pression for B derived from linear theory overestimates 
the actual value. Hence, for the rest of this paper, the 
notation B\.$ (-B2.3) signifies using a value of B divided 
by 1.5 (2.3). 

Model parameter ro, defined above, sets how 
quickly the angular momentum of particles grows after 
turnaround. This parameter is difficult to constrain an- 
alytically since the halo is nonlinear after turnaround. 
However, in Paper I we showed that a nonzero w can 
be sourced by substructure. Moreover, w influences the 
density profile at small scales since it controls how the 
pericenters of shells evolve over time. 

The trajectory of a shell after turnaround contains all 
of the velocity information in the halo. The trajectory's 
evolution equation, which follows from Newton's law, is 
given by: 



^ + (2/?-l)^+/3(/3-l)A 



2M{\) B\ 2(ru+1 _ 2 
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where £ = \n(t/t ta ) and t ta is the current turnaround 
time. The initial conditions for eq. ([5]) are A(£ = 0) = 1 
and d\/d£{£ = 0) = — j8. Calculating M(l) requires 
evolving both the shell's trajectory and A4(A) before 
turnaround [18| . Because of self-similarity, the trajectory 
A(£) can either be interpreted as labeling the location of 
a particular mass shell at different times, or labeling the 
location of all mass shells at a particular time. We take 
advantage of the second interpretation in order to numer- 
ically calculate the velocity profiles. 

Inside the shell that is currently at its second apocen- 
ter, multiple shells exist at all radii. This can be seen 
from Figure 5 of @ , which plots the location of all shells 
at a particular time. Therefore the expectation value of 
a quantity h, for example the radial velocity, at radius r 
and time t is the value of h for each shell at r weighted 
by each shell's mass. We find: 



f M ta dM, 



(h(r,t)) 



Mt, 



h(t,t*)6 D (\-X(0) 



f m a 4M*-S D (A - A(£)) 



(0) 



where M ta is the current turnaround mass, h(t,t*) rep- 
resents the value of h for the shell with turnaround time 
t*, dM* is the mass of the shell with turnaround time 
i*, and 5 D is the dirac delta function which picks out all 
shells at r. 

In order to numerically calculate h, we must relate 
h(t,t*) to A(£), the computed trajectory of the shell 
which turns around at t ta - Using self-similarity, we find: 



«t(M*) 



r ta (l+ ro -2ff)£« 



B — e 
t 
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(7) 
(8) 



where = t/t* and v t (v r ) is the tangential (radial) ve- 
locity. Using eqs. ©-([H]), and taking advantage of the 
delta function, the tangential (of) and radial (of.) veloc- 
ity dispersions become: 



_rl E i e {i ~ W+2 ^K 2 \dX/dZ\; 1 
a 2 r {r,t) = (v 2 r (r,t)) - (v r (r,t)) 2 

_ rl i>( 2 - 5 ^ K^A)M] I \d\fdt\7 1 
t* £\ e (2-30)fc idX/d^r 1 

rl ( Y.r^ im K^AV^-lrfAMi: 

t 2 



(9) 



E,e(2-3^|dAMi; 



(10) 



where is the zth root that satisfies A = A(£). In the 
above, we've imposed (vt{r,t)) = since our model as- 
sumes that the orbital planes of particles in a given shell 



3 



are oriented in random directions. Note that inside the 
shell that is currently at its second apocenter, interfer- 
ence between multiple shells forces (v r (r,t)) to quickly 
go to zero. 



III. ASYMPTOTIC BEHAVIOR 

Here, using techniques developed in Q , we analytically 
calculate the logarithmic slope of the tangential and ra- 
dial kinetic energy in two different radial regimes. We 
accomplish this by taking advantage of adiabatic invari- 
ance and self consistently calculating the total radial and 
tangential kinetic energy profiles of the halo. 

We start by parametrizing the halo mass, radial kinetic 
energy K r , tangential kinetic energy Kt and the variation 
of the apocenter distance r a . 



M(r,t) 


= K(t)r a 


(11) 


K r (r,t) 


= K r (t)r ar 


(12) 


K t (r,t) 


= K t (t)r at 


(13) 


r a 
r* 


= (w 


(14) 



In the above r* is the turnaround radius of a mass shell 
that turns around at i*. As was shown in Paper I, adia- 
batic invariance allows us to relate q and a to n. At late 
times, the orbital period is much smaller than the time 
scale for the mass and angular momentum to grow. In- 
tegrating Newton's equation and assuming n(t) and L(t) 
change little over an orbit, we find: 



2Gn(t) 
a — 1 



K--r a - 1 )-L\t){r-'-r-'). (15) 



The above relationship tell us how the pericenters r p 
evolve with time. Defining y = r a /r p and evaluating 
the above at r = r„, we find: 



y~ 2 - 1 



(16) 



In Paper I, by analyzing the radial action, we found 
that when j( < 1, n(t)r" +1 (t) — const. Therefore, for 
zu < 0, cq. (IT6l) implies that y will decrease over time. 
However, for zu > 0, K(t)r" +1 (i) = const and eq. (fT6|) 
imply that y will increase over time and will at one point 
violate y <C 1. Since we only consider bound orbits, the 
constraint y < 1 holds. At late times, as the angular mo- 
mentum continues to increase for w > 0, y ~ 1, orbits be- 
come approximately circular, the radial action vanishes, 
and L 2 (t) ~ n{t)rf l +1 [t) . Hence halos with zu < will 
have orbits that become more radial over time ()/ « 1) 
while halos with zu > will have orbits that become more 
circular over time {y ~ 1). The above insight leads to the 
following constraint. 



' ^{2w + £[a(l + n) - 3]} 
sdbj [«(!+») -3] 



if zu > 
if vo < 



(17) 



For the specific case, zu < 0, taking advantage of y <C 1, 
the adiabatic invariance arguments above, and eqs. ([1]) 
and (f?]), we can rewrite eqn. (|16[) in the form y(t,t*) — 
yo(t/t*) 1 , where: 



zu if a > 1, 

2zu/(a + l) ifa<l. 



(18) 



and 2/o r * is the pericenter of a mass shell at turnaround. 
The special case a = 1 will be addressed later. 

We next calculate the kinetic energy profiles. After 
a few orbits, shells oscillate at a much higher frequency 
than the growth rate of the halo. When calculating the 
internal mass profile, this allows us to weight each mass 
shell based on how much time it spends interior to a 
certain scale Q- Likewise, when calculating the total 
internal kinetic energy profile, we can weight each mass 
shell by both a time- averaged v 2 (or v 2 ) and a factor that 
accounts for how often the shell lies interior to a certain 
scale. For a derivation, please see the Appendix. 

Using eqs. (fTS"|) and (fTB"]) . the kinetic energy weight- 
ing Pi(r/r a ,y) at time t for a mass shell with apocenter 
distance r a , pericenter yr a , below r, is: 



Pi(u,y) = 



< v) 



Pi(u,y) = tt, — r {y<u<i) 



Pi{u,y) 



i(i,y) 
h{i,y) 
i(i,y) 



(u > 1) 



(19) 



where 



It(u,y) = 
I r (u,y) oc 
I{u,y) = 



1 {BrV " 



2 \t,r a J \t 

t 2 \r ta 
dv 



dv 



y v 2 f{v,y) 



f{v,y)dv 



f(v,y) 



(20) 
(21) 
(22) 



and 



f(v ^J^-O-*- 2 -!)) 1 ' 2 
\((v°-i-l)+A(y)( V - 2 -l)) 1/2 



if a > 1, 
if a < 1. 

The index i = {r, t} is used for shorthand to represent 
either the radial or tangential direction and the depen- 
dence of r a and r p on is implicit. Eq. ([2~Tj) is not an 
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equality since Gn(i) oc rf~ a /t 2 . Moreover, the propor- 
tionality constant varies for different radial regimes in the 
halo. 

Similar to the treatment in Paper I, self consistency 
demands that: 



recently turned around do not contribute to the kinetic 
energy inside r/r ta since we are probing scales below their 
pericenters. Mass shells only begin to contribute when 
the two argument of Pi are roughly equal to each other. 
This occurs around: 



r 

r ta . 



Ki(r ta ,t) 
M - dAL t 2 



-Pi 



M ta r 2 a \r a (t,U) 



(24) 



where dM* is the mass of a shell that turned around at 
t* and M ta is the current turnaround mass. The above 
is not an equality, even for the tangential kinetic energy, 
because of a proportionality constant, similar to 
that is not included. See the Appendix for details. Its nu- 
merically computed value does not affect the asymptotic 
slopes an. Noting from eq. @ that 



dM* 
din t* 



(3/3 - 2)M ta [ - 



3fi-2 



(25) 



using eq. ([T4|) and transforming integration variables to 
u = r/r a , we find: 



r 



ta Jrjrt 



du 



Pi(u,y(t,U)) (26) 



where k = (3/3 — 2)/(/3 — q). As u increases, the above 
integral sums over shells with smaller t*. Since the peri- 
center of a shell evolves with time, the second argument 
of Pi depends on u. The dependence varies with torque 
model (sign of w) ; hence we've kept the dependence on u 
implicit. Next we analyze the above for certain regimes of 
r/rta, and certain torquing models, in order to constrain 
the relationship between a, and k. 



A. r/Vta < 2/o, no < 

For to < 0, particles lose angular momentum over time. 
When probing scales r/r ta <C yo, mass shells with t* <C 
tta only contribute. As a result, y(t,t+) <C 1. Using eq. 
(1151), we then find: 



ya 



yo 



ur ts 



(27) 



u = yi = (yo (r/r ta ) ) 



s\ V(i+«) 



(28) 



Hence, we can replace the lower limit of integration in 
eq. (|26|) with j/i . We next want to calculate the behavior 
of eq. (|2"6l close to y\ in order to determine whether the 
integrand is dominated by mass shells around y\ or mass 
shells that have turned around at much earlier times. 
The first step is to calculate the behavior of Pi(u,y) for 
u w y. We find: 



P t {u,y)^(— X ' 
t l \ur ta 



u^il-y/u) 1 ' 2 
y-^' 2 ifa>l, 



y 



- 1 -"/ 2 ifa<l. 



(29) 



P r ( U ,y)cc^ — ) u^il-y/ufl 2 
t \ur ta J 

/ iy- 1 ' 2 if a > 1, 
x S -i+ a /2 if a< L 



(30) 



where It — 2(1 + m — q — f3)/(q — (3) and l r = a— 1. Given 
the above, we evaluate the indefinite integral in eq. (|26l) , 
noting that y is a function of u (eq. [27]) . For it ~ yi, we 
find: 



i 2 /" 



Kyi - 1) 3/2 (^) 



t 2 f du p ( ( r 
r ta J ul+k r \ ' y ° \ur ta 




l/2-k-h-a/2 



if a > 1, 
if a < 1. 
(31) 



x (u/ yi - l) 5 / 2 — 
r ta 



( l-k-l r 

' ' x J Vl 

X \ 1 2-k-lr+OL 2 -r . , 

I ' if a < 1 



if a > 1, 
1. 

(32) 



where S = l/(q—(3) and u = r/r a . For bound mass shells, 
q — (i < 0. Therefore, since 5 > 0, the first argument of Pj 
in eq. (|26[) increases while the second decreases as we sum 
over shells that have turned around at earlier and earlier 
times (u — > 00). For r/r ta yo, mass shells which most 



Following the logic in Paper I, if we keep u/y% fixed 
and the integrand blows up as yi — > 0, then the left hand 
side of eq. (|26p must diverge in the same way as the right 
hand side shown in eqs. (|3"Tj) and (|32p . Therefore, using 
eq. 
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B. 



r/r ta < yo, > 



a t - k-l t 



-S(l + k + l t )/{l + S) if col, 

-<S(l/2 + fc + Z t + a/2)(l + <5) ifa<l. 

(33) 

'5(1 -fc-Z r )/(l + 5) ifa>l, 
5(1/2 -fc-Z r + a/2) (1 + 5) ifa<l. 

(34) 



Otherwise, if the integrand converges, then the right 
hand side is proportional to (r/r ta ) li ■ Therefore, the left 
hand side must also have the same scaling, which implies 
ctj = k + Solving the above system of equations for oti 
simplifies dramatically since we have already solved for 
{a, k, q} in Paper I. Rewritten below for convenience, we 
found: 



For w > 0, the angular momentum of particles increase 
with time. As mentioned above, when probing scales 
r/r% a 2/0i mass shells with <C t ta only contribute. 
As a result, y(t,t*) ~ 1. In other words, the orbits are 
roughly circular. We can therefore replace the lower limit 
of integration in eq. (|26[) with 1 since mass shells will only 
start contributing to the sum when u ~ y ~ 1. Hence, 
the right hand side of eq. ([2T)f is proportional to {r/rta) , 
which implies on = k + Zj. Using the results from Paper 
1, reproduced below for convenience, 



= k = 



1 + n — 3nw 



, q = 2m , for < n < 3, 



(37) 



we find: 



For n < 2 



Q 



A = 



1 + n - + n) 2 + 9nm(nm - 2) 



invo 



1 + n + 3?icc7 — y+l + n) 2 + 9n7Z7(rtci7 — 2) 
ncc(4 + n) 



1 + n — 3nzu — \/(l + «) 2 + 9nuj(nzu — 2) 



3;j 



For n > 2 : 



1+n 



, q = 0. 



(35) 



Using eq. (|35[) to solve for on and making sure the solution 
is consistent, (ie: using eqs. (f33|) and (p4|) only if the 
integrand diverges as y\ — s> 0), we find: 

For n < 2 : 



12 - 9ncc7 - 2^/(1 + n) 2 + 9nvu(nzu - 2) 



1 + n + + nf + 9nw(nw - 2) 
For n>2: 

(A + n)Cinw + 2n - 10) 



a t = 



2(1 + n)(3nnj - n - 4) 
5 — n 



Q',. 



1 + n 



(36) 



The above solutions are continuous at n = 2. Tak- 
ing the no-torque limit (vo — > 0), we find at = a r — 
(5 — n)/(l + n) for all n. Assuming virial equilibrium, 
one would predict on — 2a — 1 = (5 — ri)/(l + n). Hence, 
in the no-torque limit, both the radial and tangential ki- 
netic energy are virialized. However, when w < 0, only 
the radial kinetic energy for n > 2 is virialized. All other 
profiles are out of virial equilibrium because they are 
dominated by shells which recently turned around and 
hence have not had time to virialize. Since all collapsed 
objects today have n < 2, this model predicts unvirial- 
ized halos when particles lose angular momentum after 
turnaround. 



5-n + 3nw 

a t =a r = — , for < n < 3. (38) 

1 + n — Snvo 

The no-torque case, w = 0, is consistent with the anal- 
ysis in the prior subsection. The singularity w = (1 + 
n)/3n, as discussed in Paper I, corresponds to orbits that 
are not bound. Hence we only consider w < (1 + n)/3n. 
Eq. (|3"8|) shows that the halo, for vo > 0, is in virial equi- 
librium (a; = 2a — 1). This is expected since the velocity 
profiles are dominated by mass shells that have turned 
around at t <^ t ta - 



C. y < r/rta < 1 

In this regime, we are probing scales larger than the 
pericenters of the most recently turned around mass 
shells. As a result, Pi(u,y) is dominated by the con- 
tribution from the integrand when a» y. Therefore: 



P t (u,y)<xj§- 



rr n / r 



ur ta 



if a > 1, 



-(a+l)/2 ifa< l. 



(39) 



P r (u,y) oc 



t \ur ta 



t if a > 1, 

t (a+l)/2 ifa<1 . 



(40) 



Plugging in the above into eq. ([2l)|). using the results of 
Paper I shown below for convenience, 



6 rc-2 
a = 1 , k = — — , q = — , for n < 2 



a — k 



1 + n 



4 + n 



, 9 = 0, 



3n 



for n > 2 



(41) 
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and utilizing the same divergence and convergence argu- 
ments above, we find: 



For n < 2 : 
a t = 0, 
For n > 2 : 

J° 



= 1 



-n—3nvj 



H 
n 



S < < 1±S 



if 



-'!// 



3fi 



t 



(42) 



The above solutions are continuous at n = 2. The up- 
per limit (1 + n)/3n on vd enforces that orbits are bound 
[v3 < (3/2). For n < 2, w > (5 — n)/3n results in unbound 
orbits and hence is not considered. The radial kinetic en- 
ergy follows the same profile expected from virial equilib- 
rium (a r — 2a — 1 ), even though recently turned around 
mass shells dominate the kinetic energy for n < 2. We 
believe this is a result of angular momentum not playing 
a dynamical role at these scales. Using this logic, and 
as eq. (|4"2l reveals, this is consistent with the the tan- 
gential kinetic energy not being virialized (a t ^ 2a — 1). 
Taking the limit of eq. Q36p as w — > — oo, we recover the 
same expressions as eq. (|42|) for a r . This is expected 
since in this limit, particles lose their angular momen- 
tum instantly, resulting in purely radial orbits. We do 
not expect to recover the same expressions for at since 
the tangential kinetic energy vanishes in this limit. 

Eqs. (|35l) . (|37p and (|4"Tj) determine how quickly the 
mass inside a fixed radius grows as a function of time 
(n(t) defined in eq. QT]). When q = 0, apocenters of 
mass shells settle down to a constant fraction of their 
turnaround radii, leading to a constant mass inside a 
fixed radius. For cosmologically relevant structures (n < 
1 .4), this occurs on small scales for w = 0. When q < 0, 
inward migration leads to an increasing mass inside a 
fixed radius. When q > 0, outward migration leads to a 
decreasing mass inside a fixed radius. 

This section assumed a / 1 and yet, for certain parts 
of parameter space, eqs. (j3"5"1) . (|3"T|) and (f4*Tj) give a = 1 . 
However, since the solutions are continuous as a —> 1 
from the left and right, then the results hold for a = 1 
as well. 



hat mass which characterizes the halo when it is linear 
relates to the virial mass which characterizes the halo 
when it is nonlinear. For halos today with galactic size 
virial masses, we assume the model parameter n which 
characterizes the initial density field, is set by a spherical 
top hat mass of f0 12 M Q . Specifying the top hat mass 
also sets model parameters B and p. For explicit expres- 
sions used to calculate model parameters n, p, and B, 
please see Paper I. 

When analyzing the influence of tu, we use model pa- 
rameters: n — 0.77, p = 2n, and -B1.5 = 0.39. When 
comparing to N-body simulations, we use model param- 
eters n = 0.77, p — 2n, w — 0.12, and B 15 (B2.3) = 
0.39 (0.26). This value of w ensures p oc r^ 1 on small 
scales and, as shown in Paper I, this range in B give s 
good agreement with the Einasto and NFW profiles [26| . 

For the N-body comparisons, we average p, of, and 
of in 50 spherical shells equally spaced in log 10 r over 
the range 1.5 x 10~ 4 < r/r v < 3, where r v satisfies 
M(r v ,t) = 800Trr 3 v p B (t)/3. This is the same procedure 
followed with the recent Aquarius simulation [1]. We also 
calculate r_2, the radius where r 2 p reaches a maximum. 
This radius, as well as the virial radius r Vl are commonly 
referred to in simulation papers. As discussed in Paper I, 
the density profile is isothermal for our halo over a range 
of r. Moreover, the maximum peaks associated with the 
caustics are unphysical. So, we choose a value of r_2 
in the isothermal regime that gives good agreement with 
empirical density profiles. Changing r_2 does not change 
our interpretation of the results. We find r^/fta = 0.07 
(0.05) for Si. 5 (B2.3). For reference, we find the dimen- 
sionless radius of first pericenter passage (j/q) to be 0.042 
(0.026) for B U6 (B 2 . 3 ). 

As mentioned previously, N-body simulations have fi- 
nite dynamic range. The innermost radius where the sim- 
ulation results can be trusted is set by the total number 
of particles [27] . The recent Aquarius simulations charac- 
terize their innermost radius based on the convergence of 
the circular velocity, at a particular radius, for the same 
halo simulated at different resolutions [l|. The notation 

^conv (fconv) corresponds to the smallest radius such that 
the circular velocity has converged to 10% (2.5%) or bet- 
ter at larger radii. When these radii are showed in the 
figures, we use the values quoted in Table 2 of [l| for halo 



Aq-A-2 (r& v /r_2 = 0.022 and r y c ' ' m /r- 2 = 0.052) since 
all six halos were simulated at this resolution. 



(7) 



IV. COMPARISON WITH N-BODY 
SIMULATIONS 



A. Anisotropy Profile 



In this section, using the analytic results derived above 
to gain intuition, we first analyze how w influences the 
anisotropy and pseudo-phase-space density profiles and 
then compare our numerically computed profiles to re- 
cent N-body simulations of galactic size halos [1]. 

As described in Paper I, the mass of a halo is not well 
defined when our model is applied to cosmological struc- 
ture formation since it is unclear how the spherical top 



Here we analyze the velocity anisotropy j3 v = 1 — 
of /2of for galactic size halos, where the tangential and 
radial velocity dispersions are defined in eqs. (JSJ and (JTUJ) 
respectively. Based on the analysis in Section Hill we ex- 
pect j3 v to asymptote to a constant for r/r ta <S yo since 
of oc of and (3 V to increase for yo -c r/r ta -C r v /r ta 
since of /of a r" 1 . Moreover, for radii larger than the 
first shell crossing (r ~ r v ), of = since only one shell 
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FIG. 1. The top panel shows the velocity anisotropy profile 
for a self-similar halo with model parameters n — 0.77, p = 
2n, -B1.5 = 0.39, and varying to. Smaller w leads to halos 
with more radial orbits at a particular radius. The bottom 
panel shows the smoothed velocity anisotropy profile for a 
self-similar halo with model parameters n — 0.77, p = 2n, 
vo = 0.12, and Bi. 5 (B2.3) = 0.39 (0.26). Smaller B leads 
to a larger peak width and more radial orbits. The profile is 
qualitatively similar to results from N-body simulations. The 
dimensionless radius of first pericenter passage (yo) and the 
virial radius (r v ) are labeled for clarity. The convergence radii 
for the Aquarius halo Aq-A-2 [l[ are labeled for reference. 

contributes to the dispersion. Hence, in this radial range, 
fi v = -00. 

In the top panel of Figure Q] we plot the velocity 
anisotropy for galactic size halos with varying zu. In the 
bottom panel, we plot the smoothed velocity anisotropy 
for model parameters that give good agreement with 
density profiles from simulated galactic size halos. The 
downward spikes in both panels are caustics which ex- 
ist because of unphysical radially cold initial conditions. 



In both panels, as analytically predicted, the velocity 
anisotropy asymptotes at small radii, increases at inter- 
mediate radii, and then drops off near the virial radius. 

The top panel shows that zu affects the radius of first 
pericenter passage (yo), the amplitude of f} v close to the 
virial radius, as well as the asymptotic value of f3 v at 
small radii. This behavior is intuitive since smaller values 
of zu give rise to halos populated with less circularized or- 
bits at a given radius. Note however that the envelope of 
the anisotropy profile begins to increase and become more 
radially dominated for r/r ta < yo, contradicting our an- 
alytic analysis. More specifically, for zu = —0.1, (3 V ~ 0.2 
for r/r ta < 0.001, and starts to increase at r/r ta ~ 0.001 
when it should start increasing at r/r ta ~ yo, according 
to Section IIIH This is a result of assumptions used to 
calculate a r breaking down. This is more apparent for 
zu > since the orbital period is longer. However, as 
r — > 0, the assumptions become more valid. 

In the bottom panel, we show how model parameter 
B affects the velocity anisotropy. As discussed in Paper 
I, smaller B leads to orbits that take longer to circular- 
ize and density profiles with a larger isothermal region 
(smaller yo). The bottom panel should be compared to 
Figures 9 and 10 of [I]. Though our model cannot ad- 
dress structure outside r v , the graphs are qualitatively 
very similar. The width of the peak predicted in our 
model agrees with results from N-body simulations. This 
should be expected since the parameter B was chosen so 
that the width of the isothermal region in the density 
profiles agree. However, our model over predicts the ve- 
locity anisotropy close to r v and under predicts the ve- 
locity anisotropy at small radii. In other words, at large 
scales the halo is populated with too many radial orbits 
while on small scales the halo is populated with too many 
circular orbits. 

This trend is most clearly seen in Figure [5] There we 
plot the local velocity anisotropy versus the logarithmic 
slope of the density profile for a galactic size halo with 
B2.3 = 0.26 as well as a universal relationship relating 
these two quantities that was derived by Hansen & Moore 
[2c| . The open circles correspond to 1.5 x 10~ 4 r v < 
r < r_2 while the filled circles correspond to r_2 < r < 
r v . This figure should be compared to Figure 11 of [lj. 
In the Aquarius simulation paper, the Hansen & Moore 
prediction agrees well with N-body results for r < r_2- 
However, in our Figure while there is a clear trend 
between the local velocity anisotropy and the logarithmic 
slope of the density profile, that trend does not match the 
derived relationship. Note though that our model, just 
as the Aquarius simulation claimed, does show deviations 
from the Hansen & Moore trend for r_2 < r < r v . In 
our model, this deviation is caused by a vanishing radial 
velocity dispersion. For simulated halos, other effects like 
non-sphericity or non-self-similarity may also play a role. 

The self-similar model's inability to match the ampli- 
tude of the velocity anisotropy seen in N-body simula- 
tions reveals a weakness in the model. Clearly, it is un- 
physical for all particles in a particular shell to have the 
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FIG. 2. The local logarithmic slope of the density profile 
plotted against the velocity anisotropy. The relationship re- 
lating these two quantities that was proposed by Hansen & 
Moore (2006) is also showed. Open circles correspond to 
1.5 x 10 -4 r v < r < r_2 while closed circles correspond to 
r_2 < r < r2oo- Unlike N-body simulations, our self-similar 
model does not fit the trend proposed by Hansen & Moore 
for r < r_2. This reveals a shortcoming of the model. 

same amplitude of angular momentum and the same ra- 
dial velocity. In reality, a given shell should have a radial 
velocity dispersion and should have a distribution of an- 
gular momentum that evolves with time. This possibility 
will be discussed again in Section [V] 

B. Pseudo-Phase-Space Density Profile 

Here we analyze the pseudo-phase-space density pro- 
files p/a 3 and p/a 3 f° r galactic size halos, where a 2 = 
a 2 + of. Taylor and Navarro claimed that the pseudo- 
phase-space density roughly follows the power law r^ 1 - 875 
for all halos [20j . Surprisingly, this power law matches 
predictions made by Bertschinger for purely radial self- 
similar collapse onto a spherical top hat perturbation Q . 
Ta ylo r and Navarro's claim has been verified numerically 
[U I29l - l33j . however recently the highest resolution simu- 
lations have seen evidence for departures from this power 
law near their innermost resolved radii @. 

Based on the analysis in Section IIII1 we expect the 
power law exponent to depend on {n, zu} for r/r ta <C yo- 
With model parameters {n, vu} which give p oc r _1 for 
galactic size halos, the extended secondary infall model 
predicts p/a 3 oc p/a 3 oc r -5 / 2 . This is expected for a 
virialized halo (to > 0) with p oc r _1 . For uq -C r/r ta -C 
r/r v , the model predicts p/a 3 oc r~ 2 if the radial velocity 
dispersion dominates and p/a 3 oc r -1 / 2 if the tangential 
velocity dispersion dominates. 

In the left panels of Figure [31 we plot p/a 3 and p/a 3 



for galactic size halos with varying vu. In the right pan- 
els, we plot the smoothed pseudo-phase-space densities, 
with the radius scaled by r_2, for model parameters that 
give good agreement with density profiles from simulated 
galactic size halos. In addition, we overlay the radial top 
hat solution. Scaling the radius to r_2 causes the first 
pericenter (jjo) of both models to roughly overlap, leading 
to less difference in the amplitude of the pseudo-phase- 
space density at small radii. 

The left panels show that the asymptotic slopes vary 
with vu. The numerically computed slopes match ana- 
lytic predictions. The panels for p/a 3 blow up at radii 
close to the virial radius since ay vanishes. The right pan- 
els should be compared to Figure 13 of [lj. The extended 
secondary infall model predicts that simulations of galac- 
tic size halos should see significant deviations from Taylor 
and Navarro's claim at r/r_2 ~3x 10~ 2 when analyzing 
p/cr 3 and r/r_2 ~ 10 _1 when analyzing pjo 3 . Looking 
at the residuals in Figure 13 of [lj, this prediction seems 
plausible. If higher resolution simulations do not show 
deviations from Taylor and Navarro's claim, then this 
secondary infall model would be proven incorrect since 
the model cannot consistently reproduce both the den- 
sity and velocity profiles of simulated halos. 

As shown in Section IIII1 for cosmological initial con- 
ditions (n < 2), p, a 2 and a 2 have power laws that are 
independent of initial conditions and torqueing param- 
eters in the regime yo <C r/r ta <C r v /r ta - This implies 
that the pseudo-phase-space density is universal on these 
scales. This universality on intermediate scales may have 
played a role in Taylor and Navarro's initial claim. 

Figure U] shows the difference of the pseudo-phase- 
space density power law exponent from the radial top 
hat solution, on small scales, as a function of model 
parameters n and vu. The range in n corresponds to 
10 9 < M/Mq < 10 15 . The range in vu ensures that all 
orbits are bound. According to the extended secondary 
infall model, positive vu is necessary for n > .5 in order 
to have p a r" 1 on small scales If all halos have 

p oc r^ 1 on small scales, then halos with M > 10 9 M Q 
will have p/a 3 oc r~ 5 / 2 while halos with M < 10 9 M Q 
will have pseudo-phase-space density exponents that vary 
with halo mass. If on the other hand, zu is constant for 
all halos, then as Figured] shows, the power law will vary 
with mass. 



V. DISCUSSION 

N-body simulations have revealed a wealth of informa- 
tion about the velocity profiles of dark matter halos. In 
an attempt to gain intuition for their results, we've used 
a generalized self-similar secondary infall model which 
takes into accounts tidal torques. The model assumes 
that halos self-similarly accrete radially cold mass shells. 
Moreover, each shell is composed of particles with the 
same amplitude of angular momentum. While the model 
is simplistic, it does not suffer from resolution limits and 
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FIG. 3. The left panels show p/a 3 and p/o- 3 , for a self-similar halo with model parameters n = .77, p = 2n, B1.5 = 0.39, and 
varying zu. The numerically calculated slopes match analytic predictions. First pericenter passage (t/o) for zu = and the virial 
radius (r v ) axe labeled for clarity. The right panels shows the smoothed pseudo-phase-space density profiles, with the radius 
scaled to r_2, for model parameters that give good agreement to density profiles from galactic size simulated halos. We also 
plot the radial top hat prediction. The self-similar model predicts that simulations should see deviations from the radial top 
hat power law at r/r_2 ~3x 10 -2 for p/o 3 and deviations at r/r_2 ~ 10 _1 for p/o 3 . The convergence radii for the Aquarius 
halo Aq-A-2 [fl] are labeled for reference. 



is much less computationally expensive than a full N- 
body simulation. Moreover, it is analytically tractable. 
Using this model we were able to analytically calcu- 
late the radial and tangential kinetic energy profiles for 
r/r ta < J/o and y < r/r ta <C r/r v , where y is the di- 
mensionless radius of first pericenter passage, r v is the 
virial radius, and r ta is the current turnaround radius. 

It is clear from our analysis that angular momentum 
plays a fundamental role in determining the velocity 
structure of the halo. The amplitude of angular momen- 
tum at turnaround sets the transition scale (j/o) between 
different power law behaviors in the tangential and radial 
kinetic energy profiles. Also, for collapsed objects today 



(n < 2), t*7, the parameter that quantifies how parti- 
cles are torqued after turnaround, influences the slopes 
of both the radial and kinetic velocity dispersions at small 
radii. Moreover, both the amplitude of angular momen- 
tum at turnaround and w affect the asymptotic value of 
the velocity anisotropy profile at small radii. 

For vj < 0, the self-similar halo is not virialized on 
small scales since the radial and tangential kinetic en- 
ergy is dominated by mass shells which have not had 
time to virialize. On the other hand, for w > 0, the halo 
is virialized since the dominant mass shell have had time 
to virialize. As shown in Paper I, p oc r _1 requires w > 
for M/Mq > 10 9 . Hence, positive w is favored in order 
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FIG. 4. A contour plot of dln(p/a 3 )/dlnr + 1.875, which 
shows the deviation in the pseudo-phase-space density power 
law exponent, at small scales, from the radial top hat solution. 



to reproduce N-body simulation density profiles. Quan- 
tifying zu requires analyzing N-body simulations and is 
beyond the scope of this work. However constraining w 
with simulations will provide a test for this extended sec- 
ondary infall model. 

Our model predicts that the pseudo-phase-space den- 
sity profile is universal on intermediate and large scales. 
This could potentially play a role behind the claimed uni- 
versality of the pseudo-phase-space density [20| . Since we 
do not understand how no depends on halo mass, it is im- 
possible to rule out universality on small scales, since w 
can potentially conspire to erase initial conditions. How- 
ever, if galactic size halos have p oc r -1 , then regardless 
of universality, the model predicts p/a 3 oc r -5 / 2 . While 
hints of deviations from the radial top hat solution have 
been seen in recent simulations Q , higher resolution sim- 
ulations are needed to better test the model. 

While our self-similar model has its clear advantages, 
it is also unphysical. First, all particles in a given mass 
shell have the same radial velocity. This leads to caus- 
tics. The same tidal torque mechanisms which cause a 
tangential velocity dispersion (22|, should give rise to a 
radial velocity dispersion. Second, while qualitatively 
similar, the comparison of the model's predicted velocity 
anisotropy to N-body simulation results reveals that our 
treatment of angular momentum is too simplistic. The 
model predicts too many radial orbits at large radii and 
too many circular orbits at small radii. In reality, each 
shell is composed of a distribution of angular momen- 
tum that evolves with time. In order to properly take 
these two effects into account, one would need a statis- 
tical phase space description of the halo that includes 
sources of torque. Ma & Bertschinger provided such an 



analysis in the quasilinear regime [HJ . Therefore, a natu- 
ral extension of this secondary infall model, which could 
potentially reproduce both position and velocity space 
information of N-body simulations, would be to general- 
ize Ma & Bertschinger's analysis to the nonlinear regime 
and impose self-similarity. 
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Appendix A: Deriving the Consistency Relationship 

In this Appendix, we derive eq. (|24[) . Self-similarity 
imposes that the total radial (or tangential) kinetic en- 
ergy at time t contained within radius r is given by: 

Ki(r,t) = Mta(t)^Ki(X) (Al) 

where i = {r, t} is used for shorthand to denote the radial 
or tangential direction. The kinetic energy also obeys the 
following relationship. 

i rM ta 

Ki(r,t) = - dM*v?{t,U)H[r-R(t,t*)] (A2) 
^ Jo 

where dM*,Vi(t, f*), R(t, i*) is the mass, velocity, and ra- 
dius of a shell at time t which turned around at and 
H is the heaviside function. Since, after a short time, 
shells begin to oscillate on a timescale much shorter than 
the growth of the halo, we can replace vf with a time 
averaged version (v?) and the heaviside function with 
a weighting that takes into account how often the shell 
is below r. More specifically, considering a shell with 
turnaround time such that r p (t,t^) < r < r a (t,t^), we 
have: 




where we've left the dependence on implicit. Eq. (|A3j) 
only averages over scales below r since that is where the 
shell contributes to the kinetic energy. Eq. (IA4|) is iden- 
tical to what is done in Fillmore and Goldreich, in order 
to analytically calculate the mass profile at small scales 
0. 
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Using eqs. (|A3[) and (IA4[) . generalizing to the case 
where r < r p and r > r a , plugging into eq. (|A2[) . di- 
viding by Ki{r ta ,t) and assuming a power law for the 
kinetic energy profiles in the form of eqs. (|12p and ([T^|. 



we reproduce the consistency equation. The equation has 
a proportionality constant not only because of eq. (|2"Tj) 
but also because we do not include K,i(\). This overall 
constant does not affect the asymptotic slopes. 
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